/**************************
 * 
 * Created on: 4/25/2025
 *
 **************************/

#include <math.h>
#include <stdio.h>
#include <cmath>
#include <iostream>
#include "complex.hpp"
#include "mem_ops.hpp"
#include "support.hpp"
#include "user_types.hpp"

int main(int argc, char* argv[]) {
    d_data domain_data;
    q_nums quantum_numbers;
    p_params physical_params;
    s_data solution_data;

    // Parameters
    domain_data.n_phi = 20;                    // Number of nodes in phi direction
    domain_data.n_theta = 40;                  // Number of nodes in theta direction
    domain_data.n_r = 200;                     // Number of nodes in r direction
    domain_data.R = 40.0;                      // Length of domain

    std::cout << "quantum_numbers n: ";        // Principal quantum number
    std::cin >> quantum_numbers.n;
    std::cout << "quantum_numbers l: ";        // Azimuthal number
    std::cin >> quantum_numbers.l;
    std::cout << "quantum_numbers m: ";        // Magnetic quantum number
    std::cin >> quantum_numbers.m;

    physical_params.a0 = 1.0;                  // Bohr radius
    physical_params.h = 1.0;                   // Constant
    physical_params.mp = 1.0;                  // Mass of particle
    physical_params.kp = 1.0;                  // Potential constant
    physical_params.ke = 1.0;                  // Constant in energy level equation

    // Allocate memory for solution data
    solution_data.psi = mat3D(domain_data.n_r, domain_data.n_theta, domain_data.n_phi);
    solution_data.psi_square = mat3D(domain_data.n_r, domain_data.n_theta, domain_data.n_phi);
    solution_data.r_p = new double[domain_data.n_r];
    solution_data.theta_p = new double[domain_data.n_theta];
    solution_data.phi_p = new double[domain_data.n_phi];

    // Compute wave function
    calc_psi(domain_data,
	     quantum_numbers,
	     physical_params,
	     &solution_data);

    // Print some results
    for (int i = 0; i < domain_data.n_r; ++i) {
	int theta = domain_data.n_theta / 4;
	int phi = domain_data.n_phi / 2;
	printf("real psi[%i][%i][%i]: %f, imaginery psi[%i][%i][%i]: %f\n",
		i, theta, phi, solution_data.psi[i][theta][phi].a,
		i, theta, phi, solution_data.psi[i][theta][phi].b);
    }

    printf("error in percentage : %f\n", solution_data.error);

    // Deallocate memory
    free_mat3D(solution_data.psi, domain_data.n_r, domain_data.n_theta);
    free_mat3D(solution_data.psi_square, domain_data.n_r, domain_data.n_theta);
    delete [] solution_data.r_p;
    delete [] solution_data.theta_p;
    delete [] solution_data.phi_p;

    return 0;
}
